A Simple Redistribution Vortex Method (with Accurate Body Forces) 



O.R. Tutty 

School of Engineering Sciences 
University of Southampton 
Southampton S017 1BJ, UK 



Abstract 

A circulation redistribution scheme for viscous flow is presented. Unlike other re- 
distribution methods, it operates by transferring the circulation to a set of fixed nodes 
rather than neighbouring vortex elements. A new distribution of vortex elements can 
then be constructed from the circulation on the nodes. The solution to the redistri- 
bution problem can be written explicitly as a set of algebraic formulae, producing a 
method which is simple and efficient. The scheme works with the circulation con- 
tained in the vortex elements only, and does not require overlap of vortex cores. With a 
body fitted redistribution mesh, smooth and accurate estimates of the pointwise surface 
forces can be obtained. The ability of the scheme to produce high resolution solution 
is demonstrated using a series of test problems. 
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1 Introduction 

Discrete Vortex Methods (DVM's) are Lagrangian methods for solving for rotational fluid flow in 
which the vorticity field is partitioned into a finite number of discrete vortex elements, and the 
evolution of the flow field is determined by following the motion of the vortex elements. One 
advantage of DVM's is that the computational effort is applied only in the regions of most interest. 
Also, for flow past bodies, the far field conditions are automatically satisfied as the vorticity field 
will decay to zero away from the body, avoiding the problems that can occur when truncating the 
domain in grid based methods. 

For two-dimensional inviscid flow, there is a single component of vorticity, and for inviscid flow, 
the elements are convected with constant strength at the local velocity of the fluid. For viscous 
flow, some method of modelling the viscous diffusion must to added to the numerical scheme. A 
number of schemes have been developed. One of the first applications of a DVM for viscous flow 
P] used a random walk to model the viscous effects. A random walk is simple to apply but has 
relatively low resolution and produces noisy results. A more sophisticated method is the Particle 
Strength Exchange (PSE) method, introduced in [2J, which models the viscous effects using integral 
operators. A related method is the vorticity redistribution method of Shankar and van Dommelen 
[3l [4] which, like PSE, involves redistributing circulation between the elements, but does not require 
a regular mesh, introducing elements locally as required. A more recent redistribution, based on 
time dependent Gaussian cores, is given by [5] ■ A description of many aspects of vortex methods can 



be found in the book by Cottet and Koumoutsakos 6 . Also, a comparison of four viscous methods 
which do not require the introduction of a grid can be found in [7]. 

Many papers on vortex use the total force (coefficients) to assess the accuracy of the method but do 
not present details of the pointwise surface forces. These are frequently noisy due to the irregular 
distribution of vortex particles, as can be seen in [8J, who use random walk for the diffusion. In fact, 
the total force can also be noisy and require smoothing, as in [51 [S] . In one of the benchmark papers 
using vortex methods, Koumoutsakos and Leonard |10j do present plots of the surface pressure and 
vorticity, although they observe high frequency oscillations in the surface vorticity. 

Koumoutsakos and Leonard [10] use the PSE method for diffusion. PSE involves the introduction of 
a grid to the scheme with frequent remeshing in order to maintain the accuracy of the calculation. It 
appears that the order introduced into the method by the remeshing helps smooth out the spurious 
high frequency oscillation in the forces that are found in some other methods. Hence, a possible 
way of reducing the noise would be to remesh every time step with a PSE method. However, an 
alternative is to combine the remeshing and diffusion operation in single step. This can be done 
by making the redistribution of the circulation satisfy certain constraints which have a physical 
meaning. This is a variation of the redistribution method of Shankar and van Dommelen [3] , but 
rather than redistributing circulation between vortex elements, the circulation contained in each 
discrete vortex element is independently transferred onto a small set of neighbouring nodes . The 
solution for this redistribution problem can be written explicitly as a set of algebraic equations. The 
result is a simple, efficient scheme which maintains a regular distribution of vortex elements. The 
surface force can be calculated directly from local variables, showing very good agreement with test 
data. Also, since this scheme works with the circulation of a vortex element, independent of the 
vorticity distribution in the element, there is no requirement for overlap of cores, as in the PSE 
method. 

The paper is structured as follows; first, an outline of the basic method is presented. This is 
followed by a description of the redistribution scheme. Means of calculating the body forces are 
then discussed. A brief description of the finite-spectral code used to generate test data is given. 
Test results are then presented for impulsively started flow past a cylinder and a square. Finally, 
some conclusions are given. 



2 Basic Method 

The two-dimensional incompressible Navier-Stokes equation in vorticity form is 



Dlu 8uj duj did 1 
~Dt = ~dt + ' * % ~dx + 1 'dy = Re 



" — + v iE = ^-y 2uj w 



where (it, v) are the velocity components in Cartesian coordinates {x, y), tis the time, w the vorticity 
Re the Reynolds number, and V 2 the two-dimensional Laplace operator. The velocity is normalised 
by the free stream velocity Uq, the length scales by a characteristic length L, and time by L/Uq. 
The Reynolds number is given by pLUg/fi with p and fluid density and /i its dynamic viscosity. 

Consider an individual vortex element centred on z — Zj where z = x + i y gives the coordinates in 
complex form, with a vorticity distribution 7(77) where 77 = \z — Zj\, where 2ir J Q 7(77)77^77 = 1 so 
that the distribution function has unit circulation. The vorticity field is represented by N discrete 
vortex elements so that 

- ^r,- 7 (|* (2) 



where Tj is the strength of vortex j. 
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The velocity generated by the vortex elements is given by 



N 

z — z,- 



u b +iv b = X^ r J ! _ \ 2 F{\z-z 3 \) (3) 
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where 



Ffa) = l\( S ) S d S (4) 



A number of different functions can be used for the vorticity distribution 7. In general, point 
vortices given by a delta function are not used. Instead, a smooth distribution which does not have 
the numerical problems associated with point vortices is adopted. A standard distribution, used 
here, is the Gaussian vortex, 

7(77) = — 2 e^l°* (5) 

7T<7 

where a is a measure of the core size of a vortex. This gives 

F(r,) = ^\l- e -v 2 /* 2 } (6) 

As usual, the boundary conditions at the surface of the body are satisfied by the use of a vortex 
panel method. Both standard straight, constant strength, vortex panels and the higher order panels 
given in [9] Q were investigated. The latter method uses overlapping curved panels with a linear 
distribution of vorticity along the panels. This allows an accurate representation of a smooth body 
and produces a velocity distribution which is singularity free. However, there was little difference in 
the results for the different panels for flow past a circular cylinder if enough panels were used. The 
cylinder results (Section [8]) presented below use the high order panels, but, because of the sharp 
corners, constant strength, straight panels were used for flow past a square (Section [9|. 

The velocity now consists of three components, 

u = U / +u b + u p (7) 

where u = u + i v is the fluid velocity, XJf = Uq is the free stream velocity, u;, = u b + i v b is the 
velocity generated by the vortex elements, and u p — u p + i v p generated by the panels used to satisfy 
the boundary conditions at the surface of the body. 

Numerically, an operator splitting method is used, with inviscid and viscous sub steps, satisfying 
and 

doj 1 ~ 

respectively. 

The equation for the inviscid sub step represents the fact that for two-dimensional inviscid flow 
vorticity is convected by the flow. Numerically, the element vortices are moved at the local fluid 
velocity, i.e. 

^ - u(x„t) (10) 

1 The vortex panel method presented in [5] contains a typographical error. The formula for the velocity 
generated by a panel (u*(z) in equation 16) is from a definite integral and should be the difference between 
the two terms not the sum of them. 
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A second order Runge-Kutta method is used to move the vortices at each time step 



z 



-Atufo.t,,) (11) 



z ; +1 = z» + Atu(i 3 ,t n+1/2 ) (12) 

A number of methods exist for the viscous sub step Q . The method developed here is based on the 
vorticity redistribution scheme of Shankar and van Dommelen [3J . In this method, at each time the 
circulation of vortex elements are updated through 

r n+l = ^ r n W n (13) 

k 

where W£j represents the fraction of the circulation of vortex k transferred to vortex j by diffusion 
during time step n. The summation is over a group of vortices local to Zj, the position of vortex j. 
The fractions W£j are calculated to satisfy the following constraints 

E W kj = 1 ( 14 ) 
k 

E - x ^ = E - ») = ( 15 ) 

k k 

E ( Xj - x k f = E WSjfa - Vkf = 2 hi (16) 

k k 

E W^(x 3 ~x k )(y 3 -y k ) = (17) 



where h v = i/At/Re is the characteristic diffusion distance over time At. Stability requires > 0. 
Equations ( |14|[l7 ) enforce conservation of vorticity, the centre of vorticity, and linear and angular 
momentum. 

The redistribution is performed over all vortices within a distance Rh v of vortex fc, i.e. such that 

\*j-z k \ < Rh v (18) 

The accuracy of the method depends on the value of R and a minimum value of R = 2 is required 
for a first or second order solution to exist [3] . R = ^12 was used in [3] . A solution may not always 
exist, for example if there are less than six vortices within the region (18). If no solution is found, 
new vortices are introduced a distance \/6h v from the centre (zj) until a solution exists. Further 
details of the method and its theoretical basis can be found in [3J . 



In this scheme, all vortex elements satisfying (18) must be identified. The redistribution problem 
can then be solved using a linear programming method, for example the revised simplex scheme 
found in [TTj . 



3 The Circulation Redistribution scheme 

3.1 Redistribution in Cartesian coordinates. 

The redistribution scheme above 
merits, introducing new elements as required. However, an alternative is to transfer the circulation 



T3|[T7 ) operates by transferring circulation between vortex ele- 
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onto a set of nodes at known positions. Consider a one-dimensional unsteady diffusion problem, 
with a uniform grid with grid step h. Suppose there is a vortex of strength T placed at x = x v where 
x% < x v < Xi + i, and Xj,=ih. Let 

A = {x v - Xi )/h (19) 



Then a solution of the redistribution equations is 



Si = 1-2 



h v 



- A 



(20) 



fi-x = -(1 



fi - A) 



:(1 - fi + A) 



Xk- 



where circulation T f k is transferred to the grid point x 
solution is given by 

9i+i = l-2f%V-A? 



1 



fj, 



(1 



Ai) 



9i+2 = ~(1 - ffi+i + Ai) 
with g k — for the other c/fc, and Ai = (x v — Xi+i)/h = A — 1. 



(21) 
(22) 

All other f k are zero. A second 

(23) 
(24) 
(25) 



In principle either of these solutions could be used. However, consider the case when the element is 
midway between grid points, i.e. A = ~. On physical grounds, the redistribution would be expected 
to be symmetric, but both of the solutions are asymmetric. A simple way of producing a symmetric 
solution in this case is to use the average of the two solutions, i.e. F k — h(fk+gk),k = i—l,...,i+2. 
This gives the same solution as would be obtained by using the four points i — 1 to i + 2 and assuming 
symmetry (i^_i = F i+2 and F, = F i+1 ). 



More generally, a linear combination of the two basic solutions, ( 20p2 ) and ( 23p5 ), can be used 

F k = (l-A)/ fc + Ag k: k = i-l,...,i + 2 (26) 

This produces a symmetric three point solution if the element is at a grid point, and a symmetric 
four point solution if the element is midway between grid points, with a smooth change between 
these two extreme cases. It is the simplest solution which satisfies both the redistribution equations 
and symmetry. 



Some test calcul ations have been performed for a circular cylinder case using using only three point 
formula (( 20|22 ) for x k < x a < \(x k + x k+ i) and (23 
produced undesirable short scale variations in the sol 



25) for \(x k + x k+1 ) < x < x k+1 ), but these 



ution when a vortex element moved across a 



midpoint and the redistribution changed bias. This did not occur when using the combination (26) 



Stability requires that the redistribution fractions are positive. The most restrictive case when using 



(26) is when the element is at a grid point. The condition is then 



hv 1 
~h < 72 



(27) 



As simple test problem is that of one dimensional diffusion starting from a point distribution of 
vorticity at x — xo at t = 0. Figure [l] shows the analytic and redistribution solutions for the 
vorticity for a test case at t = 1.1 with x — 1, Re = 100 and total circulation of ^(Re/n). The 
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calculation was started from the analytic solution at t = 0.1. The grid step was h = 1/100 and 
the time step was t = 0.005. Two different grids were used in the calculation. There was an offset 
between the grids and they were used alternatively at successive time steps. A number of different 
offsets, ranging from A = to A = 1) were tested, and the good agreement with the analytic 
solution shown in Figure [l] is typical. 

1.2 

1 

O.B 
3 0.6 
0.4 
0.2 



0.6 0.6 1 1.2 1.4 

Figure 1: Vorticity distribution for the one dimensional test problem at t = 1.1 with 
Re = 100, h = 1/100 and At = 0.005. The solid line is the analytical solution and the 
symbols are from the redistribution scheme. 

Consider now the two dimensional case in (x, y) with grid Xi — ih and yj — j h and a vortex element 
located at (x v ,y v ) where Xi < x v < Xi+i and yj < y v < yj+i- The one dimensional redistribution 
in the y direction is given by 

Gi = (l-S)fi + Sg h l=j-l,...,j + 2 (28) 

where 

5 = (yv - Vj)/h (29) 

and the /; and gi are obtained as in 
The two dimensional redistribution scheme is given by 

W Kl = F k G h fe = i-l,...,» + 2, l=j-l....,j + 2 (30) 

These weights satisfy all of the constraints in the original redistribution scheme 

This solution of the redistribution problem satisfies Shankar and van Dommelen condition for the 
existence of a solution; the simplest two-dimensional solution has a nine point stencil and the stability 
condition implies that the corner points of the stencil must be a distance of at least 2h v from the 
centre point. 

Figure [2] shows typical constant vorticity contours for the two dimensional diffusion problem. Again, 
two offset grids were used. The contours shown are from the numerical solution. At the scale shown 
they are identical to the contours from the analytic solution. 
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Figure 2: Vorticity distribution for the two dimensional test problem at t = 1.1 with 
Re = 100, h = h = 1/100 and At = 0.005. From the centre, the contours are for 0.8, 0.6, 
0.4, 0.2, 0.1, and 0.025. 



3.2 Redistribution in cylindrical coordinates 

The standard case of flow past a circular cylinder will be used as a test for the method. The obvious 
grid to use in this case is one using polar coordinates {r,6). However, the original redistribution 
equations must be converted to the appropriate form. Transforming between polar and Cartesian 
coordinates in standard form x = rcosO, y = rsin#, write 

x k = fa + Ar k ) cos(6» 3 + A9 k ), y k = fa + Ar k ) sin(6» i + A6 k ) (31) 



Substituting (31 1 into ( 15)17 ) and expanding, while retaining all terms 0(A 2 ) or greater, produces 
the constraints 

Wg jrj Ar k = h 2 v (32) 

fe 



kj j 



r,A0 k = 



(33) 



E 



W^Arl = J2w%r]Ael 



Y,Wl)r 3 Ar k A9 k = 



2 hi 



(34) 
(35) 



Equations ( 32p5 ) are analogous to the original constants ( l~5"fi~7 l with the Cartesian lengths replaced 
by the polar ones (Ar k and r v A9 k ) apart from the linear radial constraint (32 ) which has a nonzero 



right side. This term arises from the £ # in the polar Laplace operator; setting the right side of 



( 32 ) to zero would produce a nonphysical result in which the diffusion operator is 



c)r 2 



The redistribution formula for the radial one dimensional problem are 



Si = 1-2 



K 

Ar 



— A — 2A 



h 2 
r„Ar 



(36) 
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fi-l = 2(1 



Si + A 



fi+l — 2~(1 _ Si 



r v Ar 

hi 
r v Ar 



(37) 



(38) 



where the radial grid is given by fj 
r = r v where r 4 < r v < r i+1 . 



r + iAr, A = (r v — rj)/Ai?, and the vortex element is at 



The second redistribution solution involving the points i, i + 1 and i+2 follows in the obvious manner. 
The azimuthal redistribution solution is obtained from ( 20p4 ) with the lengths replaced by r v A9. 
Again, linear combinations of the redistribution solutions are formed, and the two dimensional 
redistribution scheme obtained as in ( 30 1 . 



3.3 Redistribution in general coordinates. 

Above the redistribution method has been presented for Cartesian and polar coordinates system. 
The method extends to other, more general, systems. For example, for a transformation given by 
x = f(rj), the condition for first moment in x becomes 

J2w^A Vj f( m ) = -tijs^j (39) 

where Ar/j = rjj — r/k and the other equations are obtained by replacing Xj — Xk by Ar/j /' (r]k)- 
Further, the redistribution problem can formulated for a general grid with fixed nodes by calculating 
the appropriate terms, a procedure equivalent to calculating the metric terms scaling the derivatives 
in the diffusion operator in a finite volume method. 



4 Computational Algorithm. 

The full scheme is a fractional step algorithm consisting of the following steps - 

1. Redistribution onto a fixed grid over a time step of \At. 

2. Convection of the vortex elements using the two-step Runge-Kutta scheme ( 11|12 ). 

3. Redistribution onto a fixed grid over a time step of \At 

A method for diffusing vorticity from the surface where it is created must be incorporated into this 
scheme. A number of methods can be found in the literature. The simplest is to create new vortex 
elements a small distance off the boundary, as in [T2|U]. In contrast, [TO] solve a diffusion problem 
using the flux of vorticity from the surface as a boundary condition. In [5] , the redistribution method 
is used to diffuse vorticity from the vortex sheet on the boundary into the interior. The last approach 
is the one used here, in the simplest manner consistent with the algorithm. The circulation from 
the vortex panels created at the start of step 2 is placed on a set of points on the boundary (the 
control points used to calculated the strengths of the vortex panels) , and these are used as new vortex 
elements added to the redistribution of step 3. All circulation distributed across the boundary during 
the redistribution steps is reflected back across the boundary, ensuring conservation of circulation 
and imposing a no flux boundary condition. 
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5 Body Forces 



The standard way to calculate the lift and drag with a DVM is to use the impulse (see e.g. [13]) - 

d f d 
(D,£) = -- / (yu,-xu)dA = -— Vl^, - Xi ) (40) 



dt 7 ^ 

i 

where T> and L are the drag and lift normalised by pU§L where L is a reference length. The 
summation is performed over the entire domain A for all circulation carrying elements. 

The surface pressure can be related to the strength of the vortex panels [HI [H] through 

Ar 

A P = (41) 
where Ar is the circulation carried by the relevant portion of the wall. 

The wall shear stress is obtained by using a finite difference formula with velocity components 
evaluated at fixed points near the surface. For the circular cylinder, the redistribution grid has 
fi = ro + (i — i) Ar where ro is the radius of the cylinder. The shear is calculated using the velocity 
at midpoints through the one sided, second order formula 

due -3iie(r ) +4ti fl (r +Ar) - u e (r Q + 2Ar) 
~dr ~ ' 2 Ar 

where ug is the azimuthal velocity. 



(42) 



6 Test Data: Flow Past a Circular Cylinder. 

There are a large number of papers which use flow past a circular cylinder as a test case. However, 
in most of these only the lift and drag (coefficients) are presented. To provide detailed data for com- 
parison, a finite difference-spectral (FDS) code was used. The streamfunction-vorticity formulation 
is used, with governing equations the vorticity transport equation (JlJ, and the Poisson equation for 
the streamfunction ip 

oj = - V 2 4> (43) 
Fourier modes were used in 9, and second order central difference formula in the radial direction. A 



one sided backwards difference formula, similar to (42), was used for the time derivative, except for 
the first time step where a backwards Euler scheme was used. The code is fully implicit, iterating 
to obtain the solution at each time step. The radial grid was stretched to give a fine grid near the 
cylinder, and place the outer boundary of the computational domain a long way from the surface. 
The non-linear terms were handled in the usual pseudo-spectral manner. 

With an impulsive start the boundary layer grows as i 1 / 2 . This scaling was used for the earlier part 
of the computation, with 

i 

r_r 0= 2 (^) V, ip = t^y, uj = f-^n (44) 



The calculation was switched to a fixed grid at t — 1, using the radial distribution from (44) at this 
time. 

This produces a relatively simple but efficient code in which high accuracy can be obtained by using 
a large number of Fourier modes and radial grid points and a small time step. Grids of up to 1024 
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complex Fourier modes, 2000 radial points, and a time step of 10 6 were used for the data presented 
below. Grid independence was checked for all Reynolds numbers. 

On the surface of the cylinder 

dp = ro_ duj_ 

d9 Re dr { ' 

which can be used to calculate the surface pressure using a reference value of zero at the front of 



the cylinder. This equation is analogous to (41 ) for the DVM. Both relate the flux of vorticity from 
the surface to the pressure gradient. 

The results produced by this code compare well with those found in other high resolution simulation 
(e.g. [10]). Also, they agree with the short time series solutions given by [T51 [TB] . 



7 Choice of grid and numerical parameters. 

As a test case for the effects of the numerical parameters and grid on the accuracy of the solution, 
the drag for the flow past an impulsively started cylinder for short time will be used. Lengths are 
scaled or the diameter of the cylinder D so that the surface of the cylinder is at r = r = 1/2, the 
vorticity is scaled by U /D where U is the free stream velocity, and the time is scaled by D/U - 
The Reynolds number is Re = Uq D jv where v is the kinematic viscosity. 

A body fitted polar grid is used in the region ro < r < r\ . This is embedded in a uniform Cartesian 
mesh for r > r\. The inner grid is arranged so that the surface of the cylinder falls midway between 
radial grid points, with = r + (i — \)Ar where Ar is the radial grid step, so that the surface 
is at r — r\/2- Azimuthally, the grid is placed at uniformly spaced points 9 = 9j with grid step 
A9 = 2ir/N where N is the number of vortex panels. The end of the vortex panels are at 9 = 9j +1 / 2: 
with the control points for the evaluation of the boundary velocity at (vq, 0j). 

The method does not explicitly allow for the t% behaviour for small time, so the solution cannot 
be expected to be accurate over the first few time steps. However, with an appropriate choice of 
parameters, solutions which achieve high accuracy after a few time steps can be obtained. 

There are six numerical parameters which must be chosen. The grid steps in r and 9 for the inner 
grid, the grid step h for the outer grid (a square grid is used here although this is not required), the 
region for the inner grid (ri), the time step At and the core size a. 

The maximum time step is fixed by the stability of the redistribution scheme. Since the algorithm 
has two redistribution substeps over time At/2, the stability condition becomes h v /h m < 1, or 
At < Reh^, where h m is the smallest of the three grid steps Ar, r$A9 and h. 

The effects of the core size were investigated in [3], and they concluded that taking a = 1/4 where 
I is the (average) length of a vortex panel was a suitable choice. This works well here also, giving 
a = r A8/4. 

The effect of the grid on the solution was investigated by fixing the number of vortex panels and 
varying Ar. Figure [3] shows the drag calculated from the impulse (40 ) for short time for an impulsive 



start with Re = 550, 400 panels, At = 0.025, and Ar = 1/200, 1/300 and 1/400. For the smallest 
value of Ar, h v /Ar w 0.85. Apart from the first few steps, the middle value (Ar = 1/300) gives 
good agreement from the drag for the FDS scheme, while for the smaller value (A = 1/400) the 
drag approaches that from the FDS scheme from below. For A = 1/200, the drag is too large. For 
clarity, only every fourth point is shown for the solutions for h = 1/300 and 1/400. All points are 
used for h = 1/200, and no filtering or smoothing has been applied for this figure. The smoothness 
is typical of the results obtained when using a body fitted redistribution mesh. 
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Figure 3: The drag for impulsively started flow past a circular cylinder at Re = 550 cylinder 
at Re = 550 with 400 vortex panels, r\ = §, At = 0.0025 and h = A: dashes, Ar = 
x, Ar = otU; +, Ar = jjL. Solid line, FDS solution. 



Figure H shows the streamwise component of the impulse, I x = J yujdA, for the same cases as in 
Figure km For an impulsive start for flow past a cylinder, the initial condition at t — + is that from 
potential flow, with a vortex sheet of (nondimensional) strength T = —2 sin(#) on the surface, giving 
I x = — tt/2. For t > 0, the impulse should decrease smoothly from the initial value. For Ar = 1/300, 
there is some (expected) irregular behaviour for the first few time steps, but the impulse is generally 
well behaved. In contrast, the smaller and larger values of Ar produce a jump in the impulse, 
followed by a relatively fast decrease for Ar = 1/200 and an increase for Ar = 1/400, consistent 
with the behaviour of the drag (Figure [3]). 

Calculations were performed for with different time steps (smaller for all three values of Ar and 
larger for Ar = 1/200 and 1/300), but the behaviour of the impulse was similar to that shown in 
Figure [4] with an overshoot for Ar = 1/200 and an overshoot for Ar = 1/300. The impulse was 
examined for a large number of other runs with Reynolds numbers varying from 150 to 9500 and a 
range of time steps, and its behaviour for short time provides a useful diagnostic as to the quality 
as the grid with regard to the ratio of the grid steps. This test could be used for other problems in 
which there is no reliable solution available to compare with, e.g. for flow past a square (Section [9] 
below) . 

For all the Reynolds numbers studied, a ratio of approximately 2/5 for the radial to azimuthal grid 
(Ar/roA#) with a = tqAO/A was found to give an accurate solution provided the time step was 
small enough. Figure [5] shows the drag for runs with At = 0.002, 200 panels and Ar = 1/150, 400 
panels and Ar = 1/300, and 600 panels and Ar = 1/450, i.e. maintaining the same scaling as for 
400 panels in Figures [3] and |4j Clearly, the grid is too coarse to provide a good match with the FDS 
solution very early in the run, but does give a reasonable value for t > 0.1. As above, there is good 
match with 400 panels, and a very close match with 600 panels. 

The outer grid step h and r\ were also varied to ensure they did not significantly affect the results 
shown in Figures [3][5j Tests were also performed with other Reynolds numbers to ensure the results 
presented below are accurate. 
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Figure 4: The streamwise component of impulse for impulsively started flow past a circular 



cylinder at Re = 550 with 400 vortex panels, At = 0.0025 and h = *, Ar = x, 
1 -■ + Ar - 1 
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Figure 5: The drag for impulsively started flow past a circular cylinder at Re = 550 cylinder 
at Re = 550 with At = 0.002; solid line, FDS solution; A, 200 panels; x, 400 panels; *, 600 
panels. 
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Figure 6: Total drag for impulsively started flow past a circular cylinder at Re = 150. Line: 



FDS solution. Symbols: DVM solution, * from the impulse(40), A from the surface forces 
(l4Tp2). 



8 DVM solutions: Flow Past a Circular Cylinder. 

Calculations were performed for impulsively started flow past a circular cylinder using Reynolds 
numbers of 150, 550, 1000, 3000, and 9500. These Reynolds numbers were chosen as they are com- 
monly used as test cases. A large of amount of test data was generated, showing excellent agreement 
with the results from the FDS code in all cases (and with data found in other studies). Represen- 
tative results are presented for three Reynolds numbers (Re = 150, 1000, and 9500), covering three 
orders of magnitude. 



8.1 Re = 150. 



Figure [6] shows the total drag obtained from the DVM and FDS methods for flow with Re = 150. 
The numerical parameters are N = 400, Ar = 1/320, At = 0.001 (h v /Ar = 0.32), n = 1 and 
h = 1/100. For the DVM, both the drag from the impulse (40) and that from the pressure and 
wall shear stress ( 4l][42 ) are shown. There is excellent agreement between all methods. Figure [7] 
shows the pressure and wall shear stress components of the drag obtained from the DVM and FDS 
schemes. Again there is excellent agreement. 

The surface distribution of pressure and wall shear stress at a single point in time (t = 1) are 
shown in FiguresJHl and [9j There is very good agreement between the values from the two numerical 
schemes. Figure [10| shows contours of the vorticity at t — 1, with the upper half of the plot showing 
the contours from the DVM method and the lower from the FDS scheme. Again, there is excellent 
agreement. 



8.2 Re = 1000. 

Figures [TT] and [12] show the total drag and drag components for the two methods for flow with 
Re = 1000. The numerical parameters are N = 600, Ar = 1/480, At = 0.0025 (hv/Ar « 0.76), 
r\ = 1 and h = 1/480. As for Re — 150, there is excellent agreement. Also, very good agreement is 



obtained for the surface forces and contours of vorticity, as can be seen for t — 3 in Figures 13 14 
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Figure 7: Drag components for impulsively started flow past a circular cylinder at Re = 150. 
Upper, pressure component. Lower, shear stress component. Lines: FDS solution. Symbols: 
DVM solution. 




Figure 8: Surface pressure p against 9 at t = 1 for Re = 150. 9 is in radians measured from 
the rear of the cylinder. The reference value is zero at the front of the cylinder {9 = ir). 
Symbols: DVM solution, line: finite difference solution. 
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Figure 9: Wall shear stress against 6 at t = 1 for i?e = 150. Symbols: DVM solution, line: 
finite difference solution. 
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Figure 10: Vorticity contours for Re = 150 at t = 1: black, positive; red, negative. Going 
from the far field towards the cylinder, the contours are for |a;| =1,2,5,10,20. The contours 
above the axis (y > 0) are from the DVM code, and those below the axis from the FDS 
code. 
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Figure 11: Total drag for impulsively started flow past a circular cylinder at Re = 10 3 . 
Line: FDS solution. Symbols: DVM solution, * from the impulse(40), A from the surface 



forces (41 42) 



and Hi 



8.3 Re = 9500. 



Flow with Re = 9500 provides a much stiffer test of the method as the flow field is more complex, 
with short scale variations in the surface forces and a much more complicated vorticity pattern than 
that found with lower Reynolds numbers. Again, however, there is very good agreement between 
the solutions for the two numerical schemes. 



Figures 16 and [17] show the total drag and drag components. The numerical parameters are N — 
3000, Ar = 1/2400, At = 0.001 {h v /Ar 0.78), r x = 6/5 and h = 1/800. The surface forces at 
t = 2 arc shown in Figures [T8[ |19| There is a high level of agreement, in particular, in the wall shear 
stress on the rear part of the cylinder where the development of relatively small scale but strong 
structures in the flow lead to large peaks and high values of the gradient along the surface. The 



complex nature of the flow can also be seen in the vorticity contours (Figure 20 1 



The values of the wall shear stress at the top and bottom shoulder of the cylinder (9 = tt/2 and 
37r/2) are slightly lower for the DVM method as compared with those for from the FDS calculations. 
However, the radial grid used for the DVM calculation near the surface is coarse as compared to that 
for FDS, and it was found that increasing the resolution gave a better match, but with an increase 
in computational effort. 



This run generated approximately 8.4 x 10 5 vortex element by t 
used by Koumoutsakos and Leonard [10]. For comparison, for Re = 
approximately 9.5 x 10 4 vortex elements at t = 3. 



= 3, a similar number to that 
150 (Figures [6p0] ) ; there were 
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Figure 12: Drag components for impulsively started flow past a circular cylinder at Re = 
10 3 . Upper, pressure component. Lower, shear stress component. Symbols: DVM solution. 
Lines: finite difference solution. 




Figure 13: Surface pressure p against 6 at t = 3 for Re = 10 3 . is in radians measured from 
the rear of the cylinder. The reference value is zero at the front of the cylinder {9 = ir). 
Symbols: DVM solution, line: finite difference solution. 
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Figure 15: Vorticity contours for Re = 10 3 at t = 3: black, positive; red, negative. Going 
from the far field towards the cylinder, the contours are for |w|=2, 5,10, 30,40. The contours 
above the axis (y > 0) are from the DVM code, and those below the axis from the FDS 
code. 
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Figure 16: Total drag for impulsively started flow past a circular cylinder at Re = 9500. 
Line: FDS solution. Symbols: DVM solution, * from the impulse(40), A from the surface 
forces (141)142). 
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Figure 17: Drag components for impulsively started flow past a circular cylinder at Re = 
9500. Symbols: DVM solution. Lines: finite difference solution. 
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Figure 18: Surface pressure p against 8 at t = 2 for Re = 9500. 8 is in radians measured 
from the rear of the cylinder. The reference value is zero at the front of the cylinder (8 = it). 
Symbols: FDS solution, line: DVM solution. 
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Figure 20: Vorticity contours for Re = 9500 
from the far field towards the cylinder, the 
contours above the axis (y > 0) are from the 
FDS code. 



at t = 2: black, positive; red, negative. Going 
contours are for |w|=5, 10,25,50, 100, 200. The 
DVM code, and those below the axis from the 



9 Flow past a square. 

There is relatively little data available for flow past a square as compared to that for a cylinder. 
However, plots of drag and vorticity contours are given in [14] for an impulsive start with Re = 100 
and the square at 15° angle of attack. A similar calculation was performed, using a uniform, body 
fitted, Cartesian grid embedded in a uniform Cartesian grid aligned the flow. Eighty one constant 
length vortex panels were used on each side of the body, and a grid step of h — 1/81 was used for 
both the inner outer grids. The change between the inner and outer grids occurred a distance 0.5 
from the body. The time step was At = 0.005, giving h v /h 0.57. With this set of parameters, the 
behaviour of the impulse early in the calculation was as expected. 

The flow at a right angled convex corner is singular, but both the pressure and the vorticity behaving 
as p-°- 456 where p is the distance from the corner [T7]. Hence, there may be large errors in the surface 
pressure obtained by integrating the panel strengths, and in calculating the lift and drag from the 
surface forces. Figure [21] shows the drag and lift obtained from both methods. There is reasonable 
agreement, given the potential for large errors. The drag is consistent with that given in [Hj. A 
calculation with 41 panels on each side of the square and h = 1/41 produced a similar result to 
that shown in Figure [21] but with a larger difference between the drag and lift calculated from the 
impulse and the surface forces. 



Vorticity contours for for t = 20 are shown in Figure 22 This figure agrees well with that in [14] , in 
particular, as regards the position and strength of the vortices downstream of the body. 

In [14] . the corners of the square were rounded to avoid unspecified numerical problems. This was 
not required for the calculations performed in the current work. 
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Figure 21: Lift and drag for impulsively started flow past a square at 15° and Re = 100. 
The top two lines are the drag, with the upper one from impulse and the lower from the 
surface forces. The bottom two lines are the lift, with the lower from the impulse and the 
upper from the surface forces. 
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Figure 22: Vorticity contours at t = 20 for impulsively started flow past a square at 15' 
and Re=100. A step of 0.5 is used with zero omitted. 
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Figure 23: Total drag for impulsively started flow past a circular cylinder at Re = 150. 
Line: FDS solution. Symbols: DVM solution using a single Cartesian redistribution grid 
with h = 0.01, At = 0.005 and 400 vortex panels. 



10 DVM scheme with a single grid. 



All of the calculations described above have used a boundary fitted redistribution mesh near the body 
and a regular Cartesian mesh further away. Care has been taken to ensure that the total circulation 
is conserved. An alternative approach, used in a number of previous studies (e.g. [TJJ 121 Q3] ) , is to 
delete any vorticity that crosses the boundary and rely on the creation process to regenerate the 
vorticity in an appropriate manner. There are several advantages to this approach. In particular, 
it allows simulation for flow past bodies of an arbitrary shape by embedding them into a regular 
grid, and simply deleting any vortex elements which are redistributed into the body. The major 
disadvantage is that the method will no longer produce high resolution values for the surface forces. 
In particular, since part of the circulation in the vortex sheet arises from the non conservative nature 



of the redistribution at the surface of the body, the surface pressure cannot be estimated using ( 41 1 
unless the deletion is accounted for. 

Simulations were performed using a circular cylinder embedded in a uniform Cartesian grid. Fol- 
lowing [5], the vortex elements created each time step were placed a distance of 1.12er above the 
surface so that the maximum velocity generated by a new element occurs at the surface, while all 
vortex elements within this distance or below the surface were deleted after the redistribution. 

The drag for a flow with Re = 150, 400 vortex panels, a time step of At = 0.005 and a redistribution 
mesh with h = 0.01 is shown in Figure [23] Also shown is the drag from the FDS scheme, showing 
good agreement with the DVM values. The vorticity distribution for both methods at t = 1 is shown 



in Figure 24 Overall there is very close agreement, although some differences can be seen near the 
surface. However, and as expected, the distribution of vortex panel strengths and the wall shear 
stress showed large high frequency oscillations. 

A further calculation was performed for Re — 150 but with 200 panels and h — 0.02 so that there 
were approximately quarter the number of vortex elements. The drag was almost the same as shown 
in Figure [23] The vorticity contours at t = 1 for both the DVM and FDS schemes are shown in 
Figure [25] Away from the body there is still good agreement between the two solutions but the lack 
of resolution near the surface with the DVM method is more apparent. 

Calculations were also performed for flow with Re = 9500. Figure [26] shows the drag for the FDS 
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Figure 24: Vorticity contours for Re = 150 at t = 1: black, positive; red, negative. Going 
from the far field towards the cylinder, the contours are for |w| =1,2,5,10,20. The contours 
above the axis (y > 0) are from the DVM code with a single Cartesian redistribution grid 
with h = 0.01, At = 0.005 and 400 vortex panels, and those below the axis from the FDS 
code. 
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Figure 25: Vorticity contours for Re = 150 at t = 1: black, positive; red, negative. Going 
from the far field towards the cylinder, the contours are for \ui\ =1,2,5,10,20. The contours 
above the axis (y > 0) are from the DVM code with a single Cartesian redistribution grid 
with h = 0.02, At = 0.02 and 200 vortex panels, and those below the axis from the FDS 
code. 
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Figure 26: Total drag for impulsively started flow past a circular cylinder at Re = 9500. 
Line: FDS solution. Symbols: DVM solutions using a single Cartesian redistribution grid: 
x, 3000 panels, h = 1/1000 and At = 0.005; *, 1600 panels, h = 1/400 and At = 0.02. 



method and the DVM scheme with two different resolutions. The better resolved DVM calculation 
has N = 3000, h = 1/1000 and At = 0.005, and the coarser calculation N = 1600, h = 1/400 and 
At = 0.02. In both cases, the redistribution grid step was chosen to approximately match the vortex 
panel length. The drag from the better resolved DVM solution shows good agreement with that 
from the FDS method, except, unsurprisingly, during the very early part of the run. The drag from 
the coarser calculation also agrees well up to t s» 1.5, but not at later times. 

Figure [27] shows vorticity contours from both the FDS method and the DVM calculation with the 
finer grid. There is good agreement, but not as close as with the body fitted grid (Figure [20| . A 
similar comparison with the lower resolution DVM solution showed a similar general structure (e.g. 
the position of the large vortices sitting off the surface) but significant differences at smaller scales, 
reflecting a lack of resolution. 



11 Conclusions 

A simple redistribution scheme for viscous flow has been presented. Unlike other redistribution 
schemes, it operates by redistributing the circulation in a vortex element to a set of fixed nodes 
rather than transferring circulation between vortex elements. A new distribution of vortex elements 
can then be constructed from the circulation on the nodes. A major advantage of the scheme is that 
the solution of the redistribution problem is given explicitly by a set of simple algebraic equations. 
A further advantage is that core overlap is not an issue for the viscous solution. 

The scheme will be stable provided the viscous diffusion length is less than the smallest mesh length 
in the problem. This restriction is similar to that found with other redistribution schemes. 

The ability of the scheme to produce high resolution solutions has been demonstrated through a 
series of test problems. Accurate estimates for both the total body and pointwise surface forces can 
be obtained when using a body fitted mesh near the surface of the body. Accurate estimates of 
the total forces on the body can be obtained when using a non conservative scheme with the body 
embedded in a Cartesian mesh. 
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Figure 27: Vorticity contours for Re = 9500 at t = 2: black, positive; red, negative. 
Going from the far field towards the cylinder, the contours are for |w|=5,10,25,50,100,200. 
The contours above the axis (y > 0) are from the DVM code with a single Cartesian 
redistribution grid with 3000 vortex panels, h = 1/1000, At = 0.005, and those below the 
axis from the FDS code. 



The solution presented for the redistribution problem is the simplest possible which satisfies the 
equations and has the required symmetry. It would be possible to obtain higher order solutions by 
extending the the computational stencil and setting higher moments to zero. There is no requirement 
to use the same mesh throughout the computational domain, and the use of local grid refinement is 
straightforward. Also, the method extends naturally to three-dimensions. 
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